//standard
#include <istream>
#include <iomanip>
#include <numeric>
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <cstdlib>

// boost
#include <boost/format.hpp>

// components
#include <tree.h>
#include <print_error.h>
#include <index_string.h>
#include <progress_indicator.h>
#include <portable_timer.h>
#include <perlutil.h>
#include <matrix.h>
#include <matrix_io.h>
#include <bcpp_parsefasta.h>
#include <stlhelper.h>
#define DEBUG_LEVEL 2
#include <debug_func.h>


#include <boost/tuple/tuple.hpp>

#include "get_arg.h"

using boost::format;
using std::vector;
using std::ostringstream;

using std::string;
using std::fill;
using std::ofstream;
using std::cout;
using std::cin;
using std::ifstream;
using boost::tie;
using std::pair;
using std::istream;
using std::max;
using std::min;
using namespace bcpp_matrix;
//using bcpp_matrix::t_matrix_f;
//using bcpp_matrix::t_dense_matrix_f;
//using bcpp_matrix::t_dense_matrix_bool;

//________________________________________________________________________________________

//	print_matrix_differences()

//________________________________________________________________________________________
void print_matrix_differences(	t_matrix_f&			tree,
								t_matrix_f&			full,
								ostream&			ostrm,
								bool				print_header,
								eCUTOFF_FORMAT		cutoff_format,
								float				cutoff,
								t_program_args&		args)
{
	vector<string> element_names = tree.names;
	t_dense_matrix_f matrix_diff_sq	  = full.matrix - tree.matrix;
	matrix_diff_sq			 *= matrix_diff_sq;

	assert(full.size() == tree.size());
	unsigned size = full.size();
	vector<unsigned> counts(size, 0);
	vector<float> diffs(size, 0.0);
	vector<float> diffs_fitch(size, 0.0);
	// M1 = tree
	if (cutoff_format == eCUTOFF_M1 || cutoff_format == eCUTOFF_BOTH)
	{
		args.VERBOSE() << "\tExclude tree relationships outside " << cutoff << " ...\n";
		for (unsigned i = 1; i < size; ++i)
			for (unsigned j = 0; j < i; ++j)
				if (tree.matrix[i][j] > cutoff)
						full.erase_element(i, j);
	}
	if (cutoff_format == eCUTOFF_M2 || cutoff_format == eCUTOFF_BOTH)
	{
		args.VERBOSE() << "\tExclude matrix relationships outside " << cutoff << " ...\n";
		for (unsigned i = 1; i < size; ++i)
			for (unsigned j = 0; j < i; ++j)
				if (full.matrix[i][j] > cutoff)
					full.erase_element(i, j);
	}


	args.VERBOSE() << "\tCalculating differences...\n";
	// calculate differences per row and fitch differences per row
	for (unsigned i = 1; i < size; ++i)
	{
		for (unsigned j = 0; j < i; ++j)
		{
			if (full.is_missing(i,  j))
				continue;
			++counts[i];
			++counts[j];
			diffs[i] += matrix_diff_sq[i][j];
			diffs[j] += matrix_diff_sq[i][j];
			if (full.matrix[i][j] > 0.0)
			{
				diffs_fitch[i] += matrix_diff_sq[i][j] / full.matrix[i][j];
				diffs_fitch[j] += matrix_diff_sq[i][j] / full.matrix[i][j];
			}
		}
	}

	args.VERBOSE() << "\tExclude tree elements without any tree or matrix relationships within " << cutoff << " ...\n";
	//
	// sort inversely by number of relationships
	// take care of ones without any relationships
	//
	std::vector<unsigned> sorted_indices;
	sort_indices(counts, sorted_indices, std::greater<unsigned>());
	order_by_sorted_indices(counts, sorted_indices);
	order_by_sorted_indices(diffs, sorted_indices);
	order_by_sorted_indices(diffs_fitch, sorted_indices);
	order_by_sorted_indices(element_names, sorted_indices);
	unsigned end = find(counts.begin(),  counts.end(),  1U) - counts.begin();
	unsigned end_1 = find(counts.begin(),  counts.end(),  0U) - counts.begin();
	for (unsigned i = end; i < size; ++i)
	{
		diffs_fitch[i]	=
		diffs[i] 		= 0.0;
	}


	args.VERBOSE() << "\t" << size << " tree elements compared\n";
	args.VERBOSE() << "\t" << end_1 << " tree elements have valid relationship pairs within " << cutoff << " ...\n";

	if (print_header)
		ostrm << "Name\tCount_pairs\tDiffs_sq_per_pair\tFitch_per_pair\n"
			     "----\t-----------\t-----------------\t--------------\n";


	//
	//	total for all pairs
	//
	args.VERBOSE() << "\tCalculating totals and averages...\n";
	float sum_squared = std::accumulate(diffs.begin(),  diffs.end(), 0.0);
	float sum_fitch_squared = std::accumulate(diffs_fitch.begin(),  diffs_fitch.end(), 0.0);
	unsigned total_count = std::accumulate(counts.begin(),  counts.end(),  0);
	ostrm	<< "Total_for_all_pairs\t"
			<< total_count << "\t"
			<< sum_squared << "\t"
			<< sum_fitch_squared << "\n";


	ostrm	<< "Avg_per_pair\t\t"
			<< std::fixed
			<< std::setprecision (5)
			<< sum_squared / total_count << "\t"
			<< sum_fitch_squared / total_count << "\n";

	//
	//	divide by counts to get averages
	//
	args.VERBOSE() << "\tNormalising difference per relationship pair count...\n";
	vector<float> counts_float(counts.begin(), counts.end());
	for (unsigned i = 0; i < end; ++i)
	{
		float inv_count = 1.0 / counts[i];
		diffs_fitch[i] *= inv_count;
		diffs[i] *= inv_count;
	}

	//
	//	Get sums
	//
	args.VERBOSE() << "\tGet total sums...\n";
	sum_squared = std::accumulate(diffs.begin(),  diffs.end(), 0.0);
	sum_fitch_squared = std::accumulate(diffs_fitch.begin(),  diffs_fitch.end(), 0.0);
	ostrm	<< "Sum_of_per_row_averages_(below)\t\t"
			<< sum_squared << "\t"
			<< sum_fitch_squared << "\n";


	if (print_header)
		ostrm << "----\t-----------\t------------\t-----------\n";



	// sort by fitch except for count = 0
	sort_indices(diffs_fitch.begin(),  diffs_fitch.begin() + end, sorted_indices, std::greater<float>());
	assert(diffs_fitch.size() == diffs.size());
	assert(diffs_fitch.size() == counts.size());
	assert(diffs_fitch.size() == element_names.size());
	assert(end == sorted_indices.size());
	args.VERBOSE() << "\tPrinting out elements with relationships within " << cutoff << " ...\n";
	for (unsigned i = 0; i < sorted_indices.size(); ++i)
	{
		assert(sorted_indices[i] < diffs.size());
		assert(sorted_indices[i] < diffs_fitch.size());
		assert(sorted_indices[i] < element_names.size());
		assert(sorted_indices[i] < counts.size());
		ostrm
				<< element_names[sorted_indices[i]] << "\t"
				<< counts[sorted_indices[i]] << "\t"
				<< std::fixed
				<< std::setprecision (5)
				<< diffs[sorted_indices[i]] << "\t"
				<< diffs_fitch[sorted_indices[i]] << "\n";
	}
	ostrm << "\n";



	args.VERBOSE() << "\tPrinting out elements with only 1 relationship within " << cutoff << " ...\n";
	for (unsigned i = end; i < end_1; ++i)
		ostrm	<< tree.names[i]	<< "\t1\t0\t0\n";
	ostrm << "\n";


	args.VERBOSE() << "\tPrinting out elements without relationships " << cutoff << " ...\n";
	for (unsigned i = end_1; i < size; ++i)
		ostrm	<< tree.names[i]	<< "\t0\t0\t0\n";


}

void LCD_matrices(t_matrix_f& full_matrix, t_matrix_f& sub_matrix, vector<string>& missing_names)
{
	vector<unsigned> sub_matrix_indices;
	vector<unsigned> full_matrix_indices;
	t_index_string indexed_names(full_matrix.size());

	for (unsigned i = 0; i < full_matrix.size(); ++i)
		indexed_names.index(full_matrix.names[i]);

	for (unsigned i = 0; i < sub_matrix.size(); ++i)
	{
		unsigned index;
		if (!indexed_names.is_indexed(sub_matrix.names[i], index))
		{
			missing_names.push_back(sub_matrix.names[i]);
			continue;
		}
		sub_matrix_indices.push_back(i);
		full_matrix_indices.push_back(index);
	}

	full_matrix.reduce_to_subset(full_matrix_indices);
	sub_matrix.reduce_to_subset(sub_matrix_indices);
}

//________________________________________________________________________________________

//	get_matrix_from_tree

//		reconstruct matrix from a tree in new hampshire format

//________________________________________________________________________________________
bool get_matrix_from_tree(const string& tree_str, t_matrix_f& tree_matrix, string& error)
{
	bcpp_tree::t_tree tree;
	bool success;
    tie(success, error) = load_tree_from_newhampshire_format_str(tree_str, tree);
	if (!success)
		return false;
	tree.to_matrix(tree_matrix);
	return true;
}

//________________________________________________________________________________________

//	read_tree_matrix_from_stream

//		reconstruct matrix from a tree in new hampshire format

//________________________________________________________________________________________
void read_tree_matrix_from_stream(istream& istrm, t_matrix_f& tree_matrix)
{
	bcpp_tree::t_tree tree;
	unsigned line_cnt = 0;
	istrm >> bcpp_tree::t_tree_io_wrapper(tree, line_cnt);
	tree.to_matrix(tree_matrix);
	tree_matrix.matrix_name = tree.get_name();
}



void read_matrix_from_stream(istream& istrm, t_matrix_f& data, t_program_args& args)
{
	t_progress_indicator dots(args.VERBOSE(), 1);
	dots.use_timer();
	unsigned cnt_duplicate_pairs = 0;
	read_matrix(istrm,  data,  args.input_format, cnt_duplicate_pairs, dots);
	if (cnt_duplicate_pairs)
		args.VERBOSE() << "\t" << cnt_duplicate_pairs << " duplicate pairs found.\n";
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main (int argc, char *argv[])
{
	t_program_args args;
	try
	{


	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	std::ios::sync_with_stdio(false);

	if (!process_command_line_options(argc, argv, args))
		return 1;


//    char buff[4096];
//    cin.rdbuf()->pubsetbuf(buff, 32768); // enable buffering

	//
	//	read matrix from input file
	//
	args.VERBOSE() << "\tReading matrix...\n";
	t_matrix_f	full_matrix;
	read_matrix_from_stream(args.matrix_strm, full_matrix, args);


	//
	//	read tree matrix from cin
	//
	args.VERBOSE() << "\tReading tree...\n";
	t_matrix_f tree_matrix;
	read_tree_matrix_from_stream(cin, tree_matrix);


	args.VERBOSE() << "\tMatching matrices...\n";
	vector<string> missing_names;
	LCD_matrices(full_matrix, tree_matrix, missing_names);


	//
	// Handle missing name from the matrix
	//
	if (missing_names.size())
	{
		args.VERBOSE() << "\tThe following names in the tree "
							"were missing from the supplied matrix:\n";
		for (unsigned i = 0; i < missing_names.size(); ++i)
			args.VERBOSE() << missing_names[i] << "\n";
		if (args.missing_matrix_as_errors)
			throw std::runtime_error("All names from the tree must be present in the matrix.");
	}


	print_matrix_differences(tree_matrix, full_matrix, cout, args.header,
								args.cutoff_format, args.cutoff, args);

	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}
	return 0;
}



